import numpy as np
import Mixing_rules as MR
from EOS_collection import CUBIC_EOS
from scipy.optimize import minimize
import time


start_time = time.time()

R = 8.314

# #输入对应态参数
# Tc = np.array([[364.8, 408.1]])
# Pc = np.array([[4.610e6, 3.648e6]])
# w = np.array([[0.148, 0.176]])
# Z = np.array([[0.275, 0.281]])
# M = np.array([[44, 44]])
# Tb = np.array([[322.02, 373.42]])

# #输入要计算的状态参数
# y1 = np.array([[0.3]])
# y2 = 1 - y1
# y = np.concatenate((y1, y2))
# x = np.ones_like(y) * 0.5
# T = np.sum(y.flatten() * Tb.flatten())
# P = 2026.5e3

#输入对应态参数 甲醇/水
Tc = np.array([[521.679, 647.1081]])
Pc = np.array([[8055865.7, 22071804]])
w = np.array([[0.56197, 0.34417]])
Z = np.array([[0.22219, 0.22968]])
M = np.array([[44, 44]])
Tb = np.array([[337, 373.42]])

#输入要计算的状态参数
y1 = np.array([[0.5]])
y2 = 1 - y1
y = np.concatenate((y1, y2))
x = np.ones_like(y) * 0.5
T = np.sum(y.flatten() * Tb.flatten())
P = 101325

EOS = CUBIC_EOS(Pc, Tc, Z, w)

def DEW(T, x):

    # 初始化
    xnew = y

    # 计算状态方程参数
    a = EOS.PR_a(T)
    b = EOS.PR_b()

    tol = 1
    while tol > 1e-5:
        # 初始化
        x = xnew
        # 计算液相逸度系数
        aL = MR.a_mix(a, x)
        bL = MR.b_mix(b, x)
        VL = 2 * bL
        VLnew = EOS.PR_V(P, VL, T, aL, bL)
        while np.sum ( np.abs(VLnew - VL) ) / np.sum ( VL ) > 1e-5:
            VL = VLnew
            VLnew = EOS.PR_V(P, VL, T, aL, bL)
        ZML = P * VL / R / T
        xa = MR.xya_mix(a, x)
        phiL = EOS.PR_PHI_MIX(P, VL, T, aL, bL, xa, b, ZML)

        # 求气相的组成
        aV = MR.a_mix(a, y)
        bV = MR.b_mix(b, y)
        VV = R * T / P
        VVnew = EOS.PR_V(P, VV, T, aV, bV)
        while np.sum ( np.abs(VVnew - VV) ) / np.sum ( VV ) > 1e-5:
            VV = VVnew
            VVnew = EOS.PR_V(P, VV, T, aV, bV)
        ZMV = P * VV / R / T
        ya = MR.xya_mix(a, y)
        phiV = EOS.PR_PHI_MIX(P, VV, T, aV, bV, ya, b, ZMV)
        
        xnew = y.flatten() * phiV / phiL
        xnew = xnew.reshape(-1, 1)
        tol = np.sum(np.abs(x - xnew))

    return (np.sum(x) - 1) ** 2

T = minimize(DEW, T, x, method = 'L-BFGS-B').x
print(f"露点温度是{T}K, {T-273.15}℃")

end_time=time.time()
elapsed_time = end_time - start_time
print(f"运行时间：{elapsed_time:.4f}秒")